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ABSTRACT 


Atmospheric turbulence has significant impairments on the operation of Free-Space Optical (FSO) communica- 
tion systems, in particular temporal and spatial intensity fluctuations at the receiving aperture resulting in power 
surges and fades, changes in angle of arrival, spatial coherence degradation, etc. The refractive index structure 
parameter C? is a statistical measure of the strength of turbulence in the atmosphere and is highly dependent 
upon vertical height. Therefore to understand atmospheric turbulence effects on vertical FSO communication 
links such as space-to-ground links, it is necessary to specify C? profiles along the atmospheric propagation path. 
To avoid the limitations on the applicability of classical approaches, propagation simulation through geometrical 
ray tracing is applied. This is achieved by considering the atmosphere along the optical propagation path as a 
spatial distribution of spherical bubbles with varying relative refractive index deviations representing turbulent 
eddies. The relative deviations of the refractive index are statistically determined from altitude-dependent and 
time varying temperature fluctuations, as measured by a microwave profiling radiometer. For each representative 
atmosphere ray paths are analyzed using geometrical optics, which is particularly advantageous in situations of 
strong turbulence where there is severe wavefront distortion and discontinuity. The refractive index structure 
parameter is then determined as a function of height and time. 
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1. INTRODUCTION 


Inhomogeneities in the refractive index of the atmosphere due to fluctuations of temperature, and to a lesser 
extent pressure, can significantly degrade the performance of free-space optical (FSO) communications links. 
Specifically the presence of turbulence can lead to spatio-temporal power fades and surges, as well as deleterious 
effects on spatial coherence, beam wander, beam broadening, etc.!~> The level of turbulence is highly dependent 
upon vertical height, therefore to properly assess vertical FSO communication links such as space-to-ground 
links, estimates of turbulence along the vertical propagation path must be specified. It has been shown previously 
that the determination of refractive index structure function C? profiles, important for the assessment of optical 
communications links, can theoretically be obtained from temperature data as measured by a microwave profiling 
radiometer, and the technique was applied to measured data with good results.° However, a complicating factor 
in application of the theoretical expressions involved is that the parameter values required in the spectral analysis 
may vary over several orders of magnitude necessitating very careful determination of the break frequency ky, 
which separates the buoyancy subrange from the inertial subrange. 


Here a computational method is applied that does not require additional parameter estimates beyond what is 
available experimentally. To accomplish this we simulate beam propagation through a turbulent medium through 
geometric ray tracing, and determine C? from cross-correlations of the ray distributions on the image plane as 
in time-lapse imagery.”*® The distance from the ray source to the image plane corresponds to the vertical height 
of interest. The propagation medium is modeled by spherical bubbles with varying refractive index deviations 
representing turbulent eddies.” Measured temperature profiles are used to generate corresponding refractive index 
distributions with which to assign values to the bubbles. In this way the fluctuations are related statistically to 
the measured profiles. 
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2. METHOD DESCRIPTION 
2.0.1 Spherical Bubble Model 


Atmospheric turbulence is a condition due to small temperature perturbations which manifest optically as fluc- 
tuations in the local index of refraction leading to intensity scintillation, and manifest mechanically as local air 
mass density fluctuations. According to Kolmogorov and Richardson, the air mass perturbations due to tem- 
perature and wind velocity are described as vorticies or eddies of varying size extending from some maximum, 
or outer scale Lo, to a minimum inner scale lg. Over this range, the inertial subrange, kinetic energy cascades 
through the turbulent eddies from Lp to Jp until finally the size scales involved are too small to have relevant 
structure and the energy is dissipated.!° 


A very effective and relatively simple method of simulating optical beam propagation through a turbulent 
medium is to model the turbulent eddies as spherical bubbles of varying refractive index®!!1° given by 


n(7s) = Nair + Ts); (1) 


where 7, is a point on or inside a sphere, n,;, ~ 1 is the index of refraction of air at atmospheric pressure, and 
(7) is a perturbation on the order of 10~4 or less. These spheres are assigned random sizes in the interval [lo, Lo] 
corresponding to the inertial subrange, and are distributed randomly throughout a volume. A distribution of 
rays, typically uniform or Gaussian, is defined at a source plane and propagated according to geometrical optics 
through the assortment of spheres until each ray intersects the image plane at some distance of interest (Fig. 1). 
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Figure 1: Bubble model.'! 


2.0.2 Geometric Ray Tracing 


The ray tracing process is straightforward, being primarily repeated applications of Snell’s law. As depicted in 
Fig. 2, each ray P is represented as a vector sum of an origin vector Pp and a unit direction vector v multiplied 
by a parameter t, oe 

P=P)+t0, (2) 
and similarly each sphere (or turbulent eddy bubble) is specified by a point on the sphere P, a center C, and a 


radius r, ar 
||P —C\|!? -r? =0. (3) 


The intersection of a ray with a sphere is determined by substituting Eq. (2) into (3) and solving for the 
parameter ¢t giving 


t=-()-C)-3 + [(-O)-o2 -(H-Or- 14, (4) 


where t must be real and the smallest non-negative value must be taken. If there are no positive real roots, then 
an intersection has not occurred. 


Figure 2: Ray-sphere intersection. 


Once an intersection has been determined, refraction through the sphere must then be considered. For 
three-dimensional applications Snell’s law is written in the vector form as 


ni (kj x fn) => nil ke x nv) (5) 


where n; and n,; are the refractive indices of the respective incident and transmitted mediums, n is the unit 
normal to the sphere at the point of intersection, # = (P — C)/||P —C'||, and k; and k; are the unit direction 
vectors of the incident and transmitted rays (Fig. 3a). 


Since Eq. (5) can be written as 


(nik; = nike) Xn= 0, 


we can define a vector F such that 


P= [rule 2) — ra(les-a)] A, 


with the property that I x f = 0, and it then follows that Eq. (5) becomes 


niki — neky = ni (ky 2) — ne (ke - fi). (6) 


From the fact that k;, - i = cos 6;, k; - fi = cos6;, and Snell’s law we have 


(a) Ray entering sphere (b) Ray exiting sphere 
Figure 3: Ray refraction through a sphere. 


which after taking the negative root (since k; - <0) and upon substitution into Eq. (6) can be solved readily 
for the refracted ray direction vector at the first interface 


hea = hy ( (") f ee A. (7) 


The refracted ray point of intersection with the back side of the sphere is 
Pp = P, = or (fa > in) kes (8) 


where P, is the first intersection point, ket is given by Eq. (7), and fi; is the unit normal at the first intersection. 
Similar to before, the ray direction vector exiting the sphere at the second interface is given by 


hea = hit (") fi fa] hm A, (9) 


where in this case n; and n; have a reversed connotation than at the first interface, ki = ket and fi is the unit 


normal at the sphere exit, # = (P, — C)/||P, — C||. The point P, becomes the new ray origin and k;,2 becomes 
the new ray direction unit vector (Fig. 3b). The process is repeated for each ray and intersecting spheres until 
all the rays have been propagated to the viewing plane. 


2.0.3 Connecting with C? 


The ray distribution on the viewing plane forms an image whose frame-to-frame variance in time can be related 
to C? through the angular tilt variance (a*), where a is the z-tilt when viewing a source in the direction 0. 
Since the frame-to-frame motion of each image is independent, the angular tilt variances (aes) in the x and y 
directions are essentially equivalent to the frame-to-frame intensity variances (I. 4) in each respective direction 
except for a scale factor related to pixel size Ap and camera focal length f. For simulation purposes this scale 
factor can be chosen arbitrarily by defining an appropriate focal length. It follows that’ 


A 
(ai, y) = 05 (Tey) (10) 


and 
(a?) = (az) + (aq) (11) 


Finally, path-averaged estimates of the refractive index structure factor C? can be determined from the relation 


L 
(a?) = of fa(z)dz, (12) 


in which L is the propagation length, and f,(z) is a path weighting function given by’ 
16 2 [oe) 2 20 1 
fa(z) = —5.82 (=) pny ier ax | | [(u cos”! u) — u? (3 -— 2u?) /(1- uw) 
at 0 o Jo 
2 
2 z ) 2 ( 24 ( z ) zd 
x c (1 z) 7 (5°) + 2u (1 tT) \ pz cos 6 


where D is the imaging aperture size, d = Ap is the pixel size, L is the path length, and z = 0 is at the image plane. 


5/6 
du dd, (13) 


3. DETERMINATION OF REFRACTIVE INDEX STATISTICS FROM 
RADIOMETER TEMPERATURE MEASUREMENTS 


Temperature and humidity profiles were taken throughout 2013 at the NASA TDRSS White Sands ground 
terminal site near Las Cruces, NM using a Radiometrics Corporation MP-3000A multifrequency microwave 
profiling radiometer (MWR) pointed to zenith (Fig. 4). This radiometer features! in total 35 calibrated channels 
of 300 MHz bandwidth each, with 21 K-band channels (22 - 30 GHz) and 14 V-band channels (51 - 59 GHz). 
Each channel has a 1.1 second integration time giving a total acquisition period of At + 40s per sample, yielding 
about 2,160 data samples per day per height. 


Figure 4: Radiometrics Corp. MP-3000A MWR at White Sands, NM. 


Since humidity fluctuations do not contribute at optical wavelengths!? only temperature profiles are consid- 
ered. Figure 5 shows an example temperature time series and Figure 6 shows a one day-averaged temperature 
profile. A linear fit gives the slope of the temperature profile as approximately -7.0 K/km. 
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Temperature Time Series for April 6th - 7th, 2013 
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Figure 5: Temperature time series at h = 1.5 km for April 6th - 7th, 2013. 
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Average Temperature Profile for 1 Day 
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Figure 6: Single day average temperature profile (April 6th, 2013). 


The fluctuating refractive index of air at position 7 and time t can be determined from temperature and 
pressure fluctuations through the relationship! 


P(r, t) 
r,t)=14+79x 10°27 14 
n(F,t) Aen (14) 
where P is the pressure in millibars and T is the temperature in kelvin. Pressure fluctuations can be considered 
negligible,'° thus P(r, t) ~ (P(r)), and pressure profiles were assumed to follow an exponential barometric model 
parameterized by available ground measurements and the temperature profile data. 


For each temperature time series at a given height, a corresponding refractive index time series was calculated 
from Eq. (14) and a 10 s moving average was also determined (Figure 7). The moving average was then removed 
(Figure 8), and a Gaussian distribution was fit to the reduced data (Figure 9). 
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Figure 7: Refractive index time series. 


Refractive Index Time Series With Moving Average Removed 
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Figure 8: Time series with moving average removed. 
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Figure 9: Histogram with Gaussian fit. 


4. APPLICATION OF THE METHOD AND RESULTS 


A single data set is comprised of about 2,160 time samples for 58 different heights from a ground height of 
1.472 km to a maximum altitude 10 km above. For each time series corresponding to a single day and height 
the refractive index distribution was determined. A source ray bundle of approximately 1 x 10° collimated rays 
with a Gaussian intensity distribution of a = 100 cm was generated and used as the source for all runs of the 
simulation program (Figure 10). The computational volume was bounded in the z and y directions to the interval 
[-5 m, 5 mJ] and in the z direction as 50 m, 100 m, or 250 m depending on the altitude interval of the radiometer. 


Figure 10: Initial ray distribution on a 512 x 512 grid. 


The source rays were propagated through 2,048 spheres with refractive indices determined by generating 
a normally distributed random number from the distribution as described above, and then adding back the 
appropriate moving average value. Example output images are shown in Figure 11. 


Figure 11: Example image plane ray distributions for h = 250 m at different times. 


In all 576 images were generated for each height. Ideally an image would have generated for each time sample, 
but this was not possible due to prohibitively long computation time. The cross-correlations between subsequent 
images was calculated from which the image shift variance was determined. After applying Eq (10) - Eq (13) with 
Ap = 0.02, D = 7.07 m, and f = 0.5 x Ap estimates of the refractive index structure factor were obtained. The 
estimates were averaged over 4 samples corresponding to the 10-min averaging of the measurement data, giving 
144 points for a particular height per day. Figure 12 shows a representative output of the algorithm at a height 
of 250 m above ground and Figure 13 shows the average C? profile for the day, plotted with the Hunfagel-Valley 
model with a value of v = 27 m/s for the high altitude wind velocity parameter.!° 
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Figure 12: C? for April 6th, 2013 at 250 m above ground level. Each point is a ten-min average. 
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Figure 13: Average C? profile for April 6th, 2013. 


The results near ground height show a clear diurnal cycle with minimum values on the order of 10~!8 m~?/3 
near sunrise and sunset, with a peak value of about 4 x 107!4 m~?/3 near midday. The profile results show a 
strong dependence of C? on altitude for heights less than 1 km, spanning nearly four orders of magnitude from 
about 107!'m~?/3 to values nearly on the order of 10~!®m~?/%, agreeing well with the Hufnagel-Valley model. 
Above 1 km the height dependence lessens before diminishing substantially above 4 km where C? values are 
generally in the range of 10716. 


5. CONCLUSIONS 


A computational method using geometrical ray tracing and techniques of time-lapse imagery was applied for 
obtaining C? profiles from microwave profiling radiometer temperature measurements. Reasonable values were 
obtained and good agreement with an established model was shown. A particular advantage of this method is 
that the specific form of C#., the temperature structure parameter, is not required - in particular for altitudes 
less than 1 km or so near the ground where the Kolmogorov “2/3 law” may not hold. Additionally, sensitive 
parameter values such as the buoyancy to inertial subrange crossover frequency are not needed. The phase-based 
nature of this computational method also has the advantage of being applicable in regimes of strong turbulence. 
Future applications of this method would benefit from computational speed-up techniques, such as implementa- 
tion on a graphics processing unit (GPU) for which ray tracing is highly suitable. 
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